arXiv: 1505.01571 vl [math.NA] 7 May 2015 


Efficient numerical calculation of drift and diffusion 
coefficients in the diffusion approximation of kinetic 

equations 

V. Bonnaillie-Noel* 1 , J.A. Carrillo^ 2 , T. Goudonffi and G.A. Pavliotis^ 2 

^Departement de Mathematiques et Applications UMR 8553, PSL, CNRS, ENS Paris 
45 rue d’Ulm F-75230 Paris cedex 05, France 
department of Mathematics, Imperial College London 
180 Queen’s Gate London SW7 2AZ, UK 
3 Inria, Sophia Antipolis Mediterranee Research Centre, Project COFFEE 
& Univ. Nice Sophia Antipolis, CNRS, Labo. J. A. Dieudonne, UMR 7351 
Parc Valrose, F-06108 Nice, France 


May 8, 2015 


Abstract 

In this paper we study the diffusion approximation of a swarming model given by a sys¬ 
tem of interacting Langevin equations with nonlinear friction. The diffusion approximation 
requires the calculation of the drift and diffusion coefficients that are given as averages of 
solutions to appropriate Poisson equations. We present a new numerical method for comput¬ 
ing these coefficients that is based on the calculation of the eigenvalues and eigenfunctions 
of a Schrodinger operator. These theoretical results are supported by numerical simulations 
showcasing the efficiency of the method. 
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1 Introduction 

In this paper we consider models from statistical physics that have the general form 


d t f + -=(«■ v r / £ 

v £ 


V/I> • W v f £ ) = -Q{f £ 

£ 


( 1 . 1 ) 


Here, f £ (t,r,v ) denotes the distribution in phase space of a certain population (of particles, indi¬ 
viduals...): r G M d , v G stand for the space and velocity variables, respectively. Equation (1.1) 


is written in dimensionless form and we consider the regime where the parameter e > 0, which 
depends on the typical length and time scales of the phenomena under consideration, is small. 
The potential (t,r) H > <3>(f, r) can be defined self-consistently, typically by a convolution with the 
macroscopic density f f £ du. On the right hand side, Q is a linear operator, which can be either 
integral or differential operator with respect to the variable v; it is intended to describe some colli¬ 
sion or friction-dissipation mechanisms. Of course, the asymptotics is driven by the properties of 
the leading operator Q: 

• We assume that Q is conservative in the sense that 

[ Q(f £ ) dv = 0 


holds. Consequently, integrating (1.1) with respect to the velocity variable, we obtain the 
following local conservation law 


d t p £ + V r • J £ = 0 


where we denote 


p £ (t,r) = j f £ (t,r,v)dv, J£ = f —f £ (t,r,v) dm 


• We also suppose that the kernel of Q is spanned by a positive normalized function v i—)■ F(v). 
As we shall detail below on specific examples, Q also satisfies some dissipation properties: 
roughly speaking these properties encode the fact that Q forces the distribution function to 
become proportional to the equilibrium F. In turn, the dissipation permits us to establish 
estimates which lead to the ansatz f £ (t,r,v ) = p £ (t,r)F(v ) + yfeg £ (t,r,v ). In order to have 
a current of order 1, we should assume J vF dv = 0, so that J £ = f g £ dv. 


Finally, let us suppose that p £ , J £ and g £ have well defined limits as e —» 0, denoted p, J, and g 
respectively. Multiplying (1.1) by y/e and letting e go to 0 yields 


Q(g) = ( v ■ V r p)F - (V r d> ■ V v F)p, 
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which provides us with a closure for the equation obtained by passing to the limit in the conservation 
law 


( 1 . 2 ) 


dtp + V r • yj vg d v J = 0. 

Owing to the linearity of Q, we rewrite this as a convection-diffusion equation: 

d t p — V r (DV r p + /CV r <hp) = 0, 
where the coefficients are defined by the following matrices 

V = - tv® Q~\vF) dv, JC = tv ® Q~\V v F) dv, 


with Q _1 the pseudo-inverse of Q, which has to be properly defined on the orthogonal of Span{F}. 
In particular, the dissipative nature of Q implies that T> is nonnegative. 

There is a huge literature on the analysis of such asymptotic problems, motivated by various 
application fields (radiative transfer theory, neutron transport, the modelling of semiconductors, 
population dynamics, etc): we refer the reader for instance to |5, [T3J |T7‘ ED, 231 EH, B7] for an 
overview of results and mathematical techniques used to handle this question, and for further 
references. The asymptotic analysis of (1.1) is of great practical interest: it is clear that the 


numerical simulation of an equation like (1.2) is by far less costly than the one of (1.1): on the 


one hand, we have eliminated the velocity variable, on the other hand, when e —> 0, (1.1) contains 
stiff terms that induce prohibitive stability conditions. However, we are left with the difficulty 
of calculating the effective diffusion and drift coefficients V and /C which relies on being able to 
calculate the inverse of the operator Q, i.e. on solving equations of the form Qcf) = h. 

We address these questions in the specific case of kinetic models for swarming for which the 
auxiliary equations that determine T> and /C do not have explicit solutions. In recent years, in¬ 
creasing efforts have been devoted to the development of mathematical models able to describe the 
self-organization of a large set of living organisms (fish, birds, bacteria...), after the pioneering work 
of Vicsek et al. pSj. Based on simple rules of information exchange between the individuals about 
their close environment, the whole population organizes itself in remarkable patterns. Starting from 
individual-based description, kinetic equations can be derived by means of mean-field regimes, as it 
is usual in statistical physics [Hi US, 155] [56]: we refer the reader to [5. [101 H2! 125, 5Tj for a thorough 
study of such models. For the models we are interested in, the interaction operator Q has the form 
Q(f) = div v (V v W(v)f + V,,/), that involves a quite complicated potential function v H > W{v). 
The equilibrium function simply reads F{v) = ^e~ w ^ v \ where Z denotes the normalization con¬ 
stant, but inverting Q is not that simple. In contrast to the standard Fokker-Planck operator 
corresponding to the Langevin dynamics with linear friction, i.e. when W(y) = v 2 /2, in general 
v i —y vF{y ) is not an eigenfunction of Q, and we do not have explicit formulas for the effective 
coefficients. Our approach for computing the coefficients is based on the spectral properties of the 
operator. In fact, a unitary transformation enables us to reformulate the Poisson equation Qcj) = h 
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in the form Hu = //, where H is a Schrodinger operator associated to a certain potential <f>, that 
depends on W; see Equation (3.1) below. In turn, the operator admits a spectral decomposition. 
Then, our method is based on expanding the functions on the eigenbasis, and then computing the 
coefficients of (p in this basis. The latter step does not present any difficulty and the computa¬ 
tional effort is concentrated on the determination of the eigenelements. In practice, we work on 
the discrete form of the equations, and we expect that only a few eigenmodes are necessary to 
capture the effective coefficients. This is the strategy we shall discuss, adopting high-order Finite 
Elements discretizations of the Schrodinger equations, which allows us to make use of performing 
computational tools [38]. 

We mention now some related problems, to which the numerical method developed in this paper 
can also be applied, and alternative approaches for the calculation of the coefficients that appear 
in macroscopic equations. First, it should be noted that a similar formalism, based on the calcula¬ 
tion of the eigenfunctions and eigenvalues of the linearized collision operator, has been developed 
for the calculation of transport coefficients using the linearized Boltzmann equation [49], Second, 
the method can be adapted to compute the effective, possibly space-dependent, coefficients that 
come from the homogenization process of advection-diffusion equations with periodic coefficients, 
see e. g. HIM], the diffusion approximation for fast/slow systems of stochastic differential equa¬ 
tions mmm, and in connection to functional central limit theorems for additive functionals of 
Markov processes [32]. In all these problems, the drift and diffusion coefficients that appear in the 
macroscopic equation can be expressed in terms of the solution of an appropriate Poisson equation. 
These effective coefficients can be alternatively calculated using either Monte Carlo simulations, 
e.g. [13 E], the heterogeneous multiscale method [IS], numerical solution of the Poisson equation 
using spectral methods, e.g. [371 EH]. Other numerical approaches include the use of linear response 
theory and of the Green-Kubo theory [30] and the expansion of the solution to the Poisson equation 
in appropriate orthonormal basis functions, e.g. Hermite polynomials. This technique, which is 
related to the continued fraction expansion [SO] has been applied to the calculation of the diffusion 
coefficient for the Langevin dynamics in a periodic potential [46] and to simple models for Brownian 
motors [33] - 

The rest of the paper is organized as follows. In Section [2] we introduce the kinetic model for 
swarming we are interested in. By using the Hilbert expansion, we detail the diffusion asymptotics, 
and we discuss some properties of the effective coefficients. The convergence can be rigorously 
justified, and a complete proof is given in the Appendix. In Section [3] we switch to the formalism 
of Schrodinger equations and we explain in detail our numerical strategy. We also provide estimates 
to justify the approximation by truncating the Fourier series. Section [4] is devoted to the numerical 
illustration of the approach, dealing with a relevant set of potentials W. In particular we bring out 
the role of the scaling coefficients that appear in the potential, with difficulties related to the so- 
called “tunnelling effect ” [EEH2E]. Section [5] is reserved for conclusions. The proof of Theorem |2.1 
on the mean held limit approximation can be found in the appendix. 
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2 Motivation: Swarming Models 


The analysis of interacting particle systems with random noise finds applications in collective be¬ 
havior and self-organization of large individuals ensembles. These models lead typically to systems 
with non-standard friction terms. In particular, the following example was proposed in the litera¬ 
ture in |15j . which includes the effect of self-propulsion and a Rayleigh type friction to impose an 
asymptotic cruising speed for individuals. This system with noise and linear Stokes friction was 
considered in m For J\f particles, it reads 


dr* = Vi dt, 


d Vi = 


(a - P \v.i\ 2 )vi 


V v ' 


X>(i r < 


D |) 


dt + V2a dTj(t), 


( 2 . 1 ) 


where T*(t) are J\f independent Brownian motions with values in and a > 0 is the noise strength. 
Here, a > 0 is the self-propulsion strength generated by the organisms, /3 is the friction coefficient, 
and 7 = a//3 is the squared asymptotic cruise speed of the individuals. Notice that (a — (3 |r>| 2 )r> = 
— aW v W(v) with W(y ) = 47 M 4 — Different models for friction can lead to asymptotic fixed 

speeds for individuals, see for instance [3B] • It is interesting, therefore, to study properties of more 
general potentials in velocity W (u), with the basic behavior of having a confinement when the speed 
is large and that zero speed is a source fixed point of the speed dynamics. The class of potentials 
for the velocity that we will consider in this paper includes, for instance, potentials of the form 
W(v) = ^|u| a - f H b , with a > b > 1. 

The mean-field limit of the stochastic particle system (2.1) above under suitable assumptions 


on the interaction potential U is given by the following kinetic Fokker-Planck equation, see [5]: 

df 


where 


+ v ■ V r f - div„ [(V r U * p)f] = di v v [a\7 v W(v)f + aW v /], 


d(^ r ) = / ’)dv, 


and * stands for the usual convolution product with respect to the space variable. 

Let us first remark that Tf — 7 is the natural relaxation time for particles to travel at asymptotic 
speed \Jot/f3. We introduce the time and length scales T and L, which are determined by the 
time/length scales of observation. They define the speed unit U = L/T that will be compared 
to V, the typical particle speed and V t h = v^T®, f^ ie typical value of fluctuations in particle 
velocity, called the thermal speed. Then we can define dimensionless variables, denoted by primed 
quantities, as 

t — T t', r — Lr', 
f{l?y,v’) = L d V d /(Tf 7 , Lr', Vv'), 
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v = Vv', 

and U'(r') = C 2 U(Lr'). 







Note that C has the dimension of velocity. With this rescaling, we obtain the following dimensionless 
kinetic equation 

df 1 

+ Vi v ■ V r f ~ V2div v [(y r U * p)f} = -Q(f), ( 2 . 2 ) 

where primes have been eliminated for notational simplicity and where the operator Q is defined 
as 


Q(f) = div„[V„W»/ + 0V„/], 


(2,3) 


with W{v) = y;H 4 — \\v \ 2 for the problem corresponding to (2.1). Here, rji, 772 , 7 , 6 and e are 


dimensionless parameters given by 


Vi = 


V 

u’ 


V2 = 


c l2 _ 
uv’ 


7 = 


a 


/3 V 2 ’ 


«=i'^ 
V 


, Tf 1 

and £= — = ——. 

T Ta 


The operator Q defined in (2.3) is the Fokker-Planck operator corresponding to the stochastic 
differential equation 

d v(t) = ~V v W(v(t)) d t + V2d dT(t), (2.4) 

where T(t) denotes standard d -dimensional Brownian motion. The generator of this diffusion 
process is 

C = -W v W{v) • V„ + 6A V . (2.5) 

Now, assume the following relation between the dimensionless parameters, with a finite asymp¬ 
totic dimensionless speed, 

T)i — 772 — e _1,/2 , 7, 6 1 ~ 0(1). 


With these assumptions, Equation (2.2) becomes 

df 


f + -= (t, ■ V,./ - div„ [(V r C/ * p )!\) = -Q(f). 


( 2 , 6 ) 


In this regime, the dominant mechanisms are the noise and the nonlinear friction. This scaling is 
the so-called diffusion scaling for kinetic equations, see 0H3] and the references therein. We remark 
that for the particular application of swarming, other distinguished limits can also be considered, 
see mi for details. Different interesting features of the model arise, depending on the relative 
magnitudes of the scaling parameters 7 , 6. 

I 11 order to obtain a closed macroscopic equation for the density p in the limit e —> 0, we use 
the standard Hilbert expansion method. Inserting the following Hilbert expansion 

f£ _ y(o) _|_ y/efW + e/( 2 ) + ... and p £ = p + a/do (1) + • • • 


into ( 2 . 2 ) and identifying terms with equal power of y/e, we get: 
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e 1 terms: Q(f^) = 0 which implies that f(°\t,r,v) = p(t,r) M(y), where 

M(v) = Z- l e- w{ ~ v)l \ Z = I e~ w{v)/0 dv, 

is the Maxwellian distribution associated to the Fokker-Planck operator Q and Z denotes the 
normalization constant. This is particularly clear by rewriting the Fokker-Planck operator Q 
as 


W) = ew L 


A/v - 1 1 


(2.7) 


e x / 2 terms: 


W (1) ) = V • V,./ (0) - div w [(V r U*p)/ ( °)] 


v ■ V r p + - V v W (■ v) ■ (V r U -k p)p 

u 


M(v). 


To invert this equation, we need to solve the following problems: 

Q(Xi) = ViM(v ) , 

1 dW 

Q( Ki) = ^WM W . 


(2.8a) 

(2.8b) 


Note that f hdv = 0 is clearly a necessary condition for the problem Q(f) = h to admit a 
solution. Assuming that the potential W increases sufficiently fast as |u| —y +cx) and given 
that the righthand side of (2.8b) is equal to — dR ^' > , from the divergence theorem we deduce 


that the solvability condition is satisfied for (2.8b). The solvability condition is satisfied 


for (2.8a) under, for example, the assumption that the velocity potential W is spherically 


symmetric. In Section 4.3, we will see how the case of nonsymmetric potentials where the 


compatibility condition (2.8a) is not fulfilled can be dealt with. 


Existence of a solution for equations of the form Q(f) = h relies on the possibility to apply 
the Fredholm alternative. For this it is sufficient to show that Q has compact resolvent in 
the space L 2 (M d ; M -1 (u)). This follows, for example, by assuming that the potential satisfies 


lim ( MV l ,IF(u)| 2 - A v W(v) ) = Too. 

|u|—»+oo \ 2 


(2,9) 


See, for example, [591 Thm. A. 19]. Under this assumption we can apply the Fredholm alterna¬ 
tive to obtain f^(t, r, v) = V r p + n- (V r U*p)p with x = (%i,..., Xd) and k = (/«i,..., Kj). 
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• e° terms: 

W (2) ) =dtf {0) + V ■ V r / (1 > - div„ [(V r U*p)fW] - div„ [(V r C/*p«)/ (0) ] , 
with p^ = f /W dv. However, using again the compatibility condition, we conclude that 

d t p — V r ■ ( W r p + IC(V r U * p)p ) = 0 , (2.10) 

where the diffusion and drift matrices are given by: 

V — — / tngixd'U and K, = — / uigiKdu. (2.11) 

,/R d 7R d 

Therefore, in the e —» 0 limit regime we expect the macroscopic density to be well approximated 


by the solution to the aggregation-diffusion equation ( |2.10[ ). We remark that by setting \ — xM 

( 2 . 12 ) 


and k = kM , the Poisson equations (|2.8|) become 

and 


£(&) = u* , 


law, , 
c(Ki) =» a :*' 0 ■ 


where £, defined by (2.5), denotes the generator of the diffusion process t (->■ v(t) defined in (2.4) 


The drift and diffusion coefficients (2.11) are also given by the formulas 

' v <g> x M dv and 


V = 


K = 


v 0 P M dv . 


This is the form of the Poisson equation that appears in functional central limit theorems and in 
the diffusion approximation for fast/slow systems of stochastic differential equations [TBs 1721 . 

Furthermore, we note that the diffusion matrix is positive def ini te and thus, we can talk properly 


about a diffusion matrix. To show this, we use (2.7) to deduce 

/ 


W)|f<to = 0 


MV, 


M{v) 


M{v) 


du. 


(2.13) 


Now, given any vector £ e \ {0}, we can compute 

£>£■£ = “/ (bM(n)]-£)(x-£) 


M(u 


du 


= -/ = 0 
./R d 


M 


X.< 

M{v) 


du > 0 


The strict equality follows from the fact that x ' £/Af ^ const, as we can immediately deduce 


from (2.8a). 


The analysis of the asymptotic regime remains technically close to the derivation of the diffusion 
regimes for the Vlasov-Poisson-Fokker-Planck equation nmm muss]. We detail in the Appendix 
the proof of the following statement. 

















Theorem 2.1 Let us consider a sequence of initial data ff nit > 0 that satisfies 

SU P // (1 + |r| + W + | ln(/f nit D)/ T £ nit dv dr = M 0 < 0. 


£>0 


We suppose that the potentials U and W satisfy the technical requirements listed in Appendix \A\ 
Let 0 < T < oo. Then up to a subsequence, still labelled by e, the associated solution f £ to ( 2.6[ ) 
converges weakly in ^((CfT) x R d x R d ) to p(t,r)M(v), p e converges to p in C°([0, T]; L^lP) — 


weak), with p being the solution to (2.10) and initial data p{t = 0,r) given by the weak limit in 
L\R d ) of f ff nit dv. 

For the time being, let us discuss the numerical evaluation of the effective coefficient matrices 


T>,/C. In the particular case considered here, the right hand side of both equations in (2.8) (or 
equivalently ( 2 . 12 | )) is of the form v times a radial function, then we can simplify the diffusion and 
drift matrices by taking into account the symmetries of the problem. We leave the reader to check 
the following result. 

Lemma 2.1 Given v G R d , let us define for any indices i, j the linear operators 

Tijiy ) : Vi Vj and Ti{v) : v { U -v^ 

that exchange V{ and v 3 and change v, by —i\ respectively. Then, the following relations hold: 


Xi{%j(v)) = Xj{v) and Xi(%{v)) = ~Xi ( v ), 


and 


Ki(Tij(v)) = Kj(v ) and Ki(Ti(v)) = -tti(v). 
As a consequence, we deduce that there exist reals D > 0 and k such that 


T>ij = - ViXjiy) dv = D6. 
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and 


ICij — 


ViKj(v) dn = kS ij , 


and the macroscopic aggregation-diffusion equation (2.10) becomes 


d t p = DA r p + kV r ■ ({S7 r U x p)p) . 


(2.14) 


Therefore, in order to compute the effective macroscopic equation (2.14), we need to find the 


solutions to one component of (2.8) (or (2.12)). They are given by explicit formulae in very few 
cases. A particular example is given by the quadratic potential W{v) = then x{ v ) = dn{v) = 
—vM{y). This is due to the fact that — vM(v ) is an eigenfunction of the Fokker-Planck operator 
Q with quadratic potential. However, for more general, non-quadratic potentials such as the one 
used in the swarming model, it is not possible to obtain explicit formulas for the coefficients of the 


limiting equation (2.14). In the next section we will study this problem by eigenfunction expansion 


of the Fokker-Planck operator and we will discuss how to accurately approximate those coefficients. 
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3 Approximation of the Diffusion and Drift Coefficients 

We start by recalling a well-known (unitary) equivalence between Fokker-Planck (in L 2 (R d ] M -1 dn)) 
and Schrodinger operators (in L 2 (R d )). By setting 



it is easy to check that TL reduces to the Schrodinger operator 

TL{u) = —OAu + $(v)u, 

with the potential 

<b(v) = - i -b v W(v) + f § \V v W(v)\\ 

We remark that <P is precisely the potential that appears in Assumption 
—9 A + $(u) is dehned on the domain 

D(H) = {ue L 2 (R d ), G L 2 (R d ), A u G L 2 (R d )}. 

We point out that this also defines the domain of the operator Q. When working with the operator 
TL, the Lebesgue space L 2 (R d ) is a natural framework; it what follows, we shall denote by (•(•) the 
standard inner product in L 2 (R d ). Using classical results for Schrodinger operators, see for instance 
[481 Theorem XIII.67], we have a spectral decomposition of the operator TL under suitable confining 
assumptions. 

Lemma 3.1 We suppose that $ G Lj oc (R d ) is bounded from below, and satisfies <P(u) —>■ +oo 
as |u| —>• oo. Then, "H -1 is a self-adjoint compact operator in L 2 (R d ) and TL admits a spectral 
decomposition: there exist a non decreasing sequence of real numbers {A n } nS N —* oo, and a L 2 (R d )- 
orthonormal basis {T n } ngN such that 1-L(f$> n ) = AA 0 = 0, Ai > A > 0. 


(3,1) 


2.9 The operator TL = 


We remark that the spectral gap A of the Schrodinger operator TL is the Poincare constant in 
the Poincare inequality associated to the Fokker-Planck operator Q, i.e. the Poincare inequality 
for the probability measure M(y) dv = Z~ 1 e~ w( ' v ^ e dv. 


Notice that the property (2.13) implies that the operator TL is positive and that the kernel is 
spanned by \[M. We wish to solve the equations 


Q( X ) = vM and Q(k) = -V v WM , 
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or equivalently, 


TL(u x ) = —v\[M with u x = 


X 


K 


, _ and TL(u k ) = — -V v W \[M with u K = ,_ 

VM 9 yfM 
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The diffusion and drift coefficients are defined by the quadratic quantities 


V = — / v 0 x dv = / - H(n x ) <8) u x dv 


(3,2) 


and 


/C = 


r;ig)fi:dn= / r H{u x )®u K dv. 


(3.3) 


The Schrodinger operator is selfadjoint in L 2 (R d ) and under Assumption (2.9) it has compact 
resolvent. We denote its eigenvalues and eigenfunctions by {A n , T n }^. A lot of information 


on the properties of the eigenvalues and eigenfunctions of Schrodinger operators is available 
Using the spectral decomposition of T~L we can obtain a formula for the effective drift and diffusion 
coefficients. These formulas are similar to the Kipnis-Varadhan formula for the diffusion coefficient 
in the functional central limit theorem for additive functionals of reversible Markov processes |31| . 
We can use these formulas to develop a numerical scheme for the approximate calculation of the 
drift and diffusion coefficients. Indeed, if we develop 'H(u x ), u x , T-L{u n ) and u K in the eigenbasis, 
we get 

OO OO 

Vk, 

A/c 


'HM = U x = (m x |To)T 0 + 


fc=i 


k=1 


and 


with 


H(u K ) = , u K = [u K |T 0 )Tq + ^ , 

k 


k=1 


k=1 


Vk — 'H(u x )d> k dv and cu k = 'H{u K )'S’k du, for all k G N. 


Substituting in (|3.2|) and (J3.3[) , we obtain the following formulas for the diffusion and drift matrices: 

(3.4) 




Vk ®Vk 


k= 1 


A i 


and 




Vk <E> Uj k 


k=1 


At 


We can obtain approximate formulas for the drift and diffusion coefficients by truncating the data: 
given h G L 2 (R d ) and denoting 

OO p 

h = ( k d>k , c k = hd> k dv , for all k G N, 

fc =0 ^ 
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We denote by u N the solution of l-L{u N ) = h N , and compare it to u, the solution of T-L(u) = h. We 
use the Sobolev and Cauchy-Schwarz inequalities as follows 

A\\u — u N \\ 2 L 2 < {H(u — u N )\u — u N ) — (h — h N \u — u N ) < ||h — h N \\ L 2 \\u — u N \\ L 2 . 


Hence, we get 

||w - u n \\ L 2 < i ||h - h N \\ L 2. 

The accuracy of the approximation is therefore driven by the accuracy of the approximation of the 
data h by its truncated Fourier series: we need information on the behavior of the eigenvalues X n 
for large n and on the accuracy of the spectral projection of h. We proceed by analogy to the 
standard theory of Fourier series, where the behavior of the Fourier coefficients is related to the 
regularity of the function. More precisely, the estimate 

\\h-h N \\ L 2<CX- N k +1 (3.5) 

holds for some k > 0. Assume 'H fc (/r) e L 2 (M d ); by the spectral decomposition, we obtain 

oo oo 

K k (h) = J2(K k (h)\y n )y n with \(U k (h)\^ n )\ 2 < oo. 

72=0 72=0 

Now, we estimate the difference as 

\\h - ft" \ih = E i('*i«'»)i 2 = E 

n>N n>N A?l 


where we used that - H( v h n ) = X n ^ n and the fact that % is self-adjoint. Therefore, we deduce that 


|| h — h 


N || 2 


< 


I L 2 — 12 k 

Ar + 1 n>N 


E \CH k (h)\*r 


< 


1 


\2k 

a N+1 




since the eigenvalues are in increasing order, leading to the desired estimate (3.5). A similar 
argument shows that if H k+1 (h) e L 2 (M d ), then 


\mh- h N )\\ L . < ^wH*+\h)\\u. 


k+li 


\k 

a N+1 


A direct application of the strategy above to 


1 . 


v h(v) = —v\JM(v) and to V ^h(v) = - -V v My/M(v), 

0 

which satisfy H k {h) G L' 2 (R d ) for all k > 0, together with the symmetry of the potential in Lemma 
(2.1) leads to the main result of this section estimating the error due to the truncation in (3.4). 
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Theorem 3.1 Given V N and K, N the truncated diffusion and drift coefficients defined by 


V N = I H(u*) ®u*dv = D n I, and IC N = j H{u*) dv = K N I, 

with 'H(u^)(v) = (—Vy/M(v)) N and / H(u^)(v) = (—j}V v M^M(v)) N , or equivalently 


X N U.I2 


N 


dN = 7,H 


\vk I' , t'N Is-^Vk-Uk 


d Afc 

K =1 


and K N = - 


d “ A k 

k= 1 


(3,6) 


then the following error estimate holds: for all k > 0 and all N e N, there exists > 0 (depending 
on k but not on N) such that 


\D-D n \ + \K-K n \ < C k A 




4 Numerical Approximation and Simulations 


The numerical method in practice works as follows: 


• Step 1.- % —* 'H R \ We consider the problem set on [— R, +R] d , with R 1 and completed 
with homogeneous Dirichlet boundary conditions: owing to the functional framework, we 
expect that the eigenfunctions are localized around the origin and are exponentially decreasing 
far away the wells of the potential, see [1]. In particular, it holds under our assumptions on 
the behavior of the potential as |u| —> +oo, from [29, Tlnn. 3.4, Thin. 3.10]. Then, we choose 
R large enough to reduce the truncation error. In our examples, we fix R = 10. 


Step 2.- H b —» We use Finite Elements methods to discretize the operator R r . In 

practice, we use the library MELINA [38], a uniform mesh of [—1?,/?] with 1000 uniform 
Pio elements and a quadrature of degree 21. We have chosen this method since we need an 
approximation of the solution of the PDE with high order accuracy, see Remark |4.2[ Here and 


in what follows we denote by p, > 0 a measure of the accuracy of the underlying discretization 
method. It thus contains the information both on the refinement of the mesh, and the degree 
of the piecewise reconstruction. 


• Step 3.- (A n , T n ) —» (A^, Having at hand the discrete operator on the truncated 

domain, denoted we determine its first N eigenclcments ((Af’ M , Tf’ M ),..., (A^’ M , T^)). 

In the Finite Elements framework, the eigenvectors d/ R,tl are piecewise polynomials functions 
approximating the nth eigenfunctions. We recall that the eigenvectors form an orthonormal 
family. 
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Step 4- ( r] n ,uj n ) —> (rj R,fl , Given the data 


T-i(u x )(v ) = -v\/M(v) and U{u K ) (v) = --V v W(v)y/W(v), 

U 

we compute the corresponding N Fourier coefficients by using an appropriate quadrature 
formula, depending on the approximation framework, for the discrete analogue (r) R ’^, <^’ M ) of 

Vn= I H(u x )(v)^ n (v)dv, u n = j H(u K ){v)4f n (v) dv. 

Our results here are computed using a simple composite rectangular rule. 

Step 5.- (D, K ) — y (D R,,J ' ,N , K R,IX ’ N ): We now approximate the diffusion and drift coefficients 


using (3.6) to conclude 

d R,^n = ± 
4 


N 


1 |^| 2 

“ n =1 


JV 


and 


cl 




Un' ■ Pn 


n= 1 


\R,fi 

An 


(4,1) 


Once the eigenelements are known, the computational cost of the evaluation of the coefficients 
{j] R)il , uj r ’^) is linear with respect to the size of the linear problem to be solved (that depends 
directly on /i). Hence, the main source of the computational cost relies on the determination of the 
N eigenpairs. 

For the potentials that we consider in this paper, Lanczos-like algorithms can be used. As an 
iterative method, its computational cost cannot be estimated a priori. Nevertheless, we expect that 
only a few eigenpairs can provide an accurate result (for the quadratic case, the problem is exactly 
solved with the first eigenpair associated with a positive eigenvalue), so that the resolution would 
be far less costly than solving the linear system, a problem that, for small /i’s, would also require 
iterative methods. Here we use standard Lanczos techniques; we refer the reader to [M, (521 El] for 
further information on these methods and to [6j for the computation of the first few eigenpairs of 
complicated Schrodinger operators, based on Finite Elements approximations. 

We will show numerical simulations for three different potentials in one dimension given by: 

• Case A.- The symmetric smooth potential given by W{v) = Ay 4 — \v 2 with 7 > 0. 

• Case B.- The symmetric singular potential given by W(y) = A-u 4 — ||n| 3 with 7 > 0. 

• Case C.- The tilted smooth potential given by W{v) = A v 4 — \v 2 — 8v with 7, 5 > 0. 

We can gather all of them in a single potential 


W(v) = _ ||„|3 


cr 


5v 


with 7 > 0, 8 > 0, crG {0,1}. 


(4.2) 


14 








Note that in Case C (or a — 0, 5 > 0 in (4.2)), the potential W is not symmetric and the 
compatibility condition for solving the auxiliary equation (2.8a) is not satisfied. We shall see how 
the theory can be adapted to this case (see Section 4.3). 


Remark 4.1 Note that in the one dimensional case, the drift coefficient can be expressed in a 
simpler form. This is due to the fact that we can solve explicitly the one dimensional Poisson 
equation, up to quadratures 


Sec. 13.6]. Indeed, according to (2.12) and the expression of the 
operator Q in (2.7), we notice that 


Id W /s „, s dM, , „ d / A {if 

Q(if) — 77 —;— ( v)M(v ) =--—(v) = 9— ( M— ( — 

9 dv J ' J dv K J dv V d v\M 


Then, direct integration yields 


V’O) = 0\ ~ v + 


>M ( v ) dvJ M[y) dv. 


Therefore the drift coefficient defined in (2.11) becomes 




-v+ vM(v)dv) M(v) dv. 


(4.3) 


In cases A and B, we have J R vM(v) dv = 0 by symmetry and K is given by K = | f R v 2 M(v) dv. 
This explicit formula will be used to check the accuracy of the method. 

In order to reduce the number of free parameters, we rescale the velocity by defining v = ^/yv 
into the Fokker-Planck operator in (2.7), to get 


Q(f) = div„ 


V v W(v)f + —V„ / 

7 


where we have dropped the tildes for notational simplicity, with the rescaled potential 




W(v) = -v 4 — 3 


v| 3 -v 2 -v with 7 > 0, S > 0, a G {0,1} 

2 V7 


In this way, 0—^0 and 7 —» 00 play the same role. I11 fact, for the symmetric smooth potential of 
Case A (a = 0), all terms involving 7 disappear in the rescaled potential and we can remove one 
parameter by setting 9—1. In our simulations, we consider the Schrodinger operator 

~ 9 d 2 u ~ . 

U{u) = —— + $(v)u, 

7 Oir 
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with the potential 


$ 0 ) = — 



( 3y2 - 2a Vl\ v \ - (! - ^)) + ^ - <Ty/lv\v\ - 


(1 - a)v - 

\/7 


in each of the different cases above and denote by Xki'y) the fc-th positive eigenvalue. 

The first eigenvalue for R is simple, equal to 0. 

In Cases A and B, the potentials are spherically symmetric. Then the second eigenvalue Ai(y) 
tends exponentially to 0 as 7 tends to 00 (ED, E3J and [531 Thm. 1.5]. In fact, a careful reading 
of these references gives that Ai(y) ~ e -C7 for some positive constant c and A 2 (y) — 0(1). This 
is a manifestation of the tunnelling effect. This behavior leads to numerical difficulties. Indeed, 
we have to capture two simple eigenvalues (0, Ai( 7 )) but with Ai(y) —* 0 exponentially fast as 
7 —> 00 . The first eigenfunction is symmetric and the second one antisymmetric. Numerically, 
when 7 is very large, the gap between 0 and Ai(y) becomes negligible compared to the order of 
the accuracy of the method or even compared to machine precision. Then numerically it appears 
as if the problem has a double eigenvalue. Then the computation breaks down. Similar difficulties 
appear for the magnetic tunnelling effect, see [7]. A good way to determine whether or not the 
computation is accurate is to look at the eigenfunction: as soon as the symmetry is broken for the 
first two eigenfunctions, the computation is wrong. 


Remark 4.2 We have also tested the method by using the standard Finite Difference discretization. 
We roughly obtain similar results for small values of 7 ’s and with the same number of numerical 
unknowns as for the Finite Elements algorithm (which means with a very refined grid for the finite 
discretization method). Discrepancies appear as 7 increases: the loss of symmetry of the eigen¬ 
functions is sensitive earlier. This is reminiscent to the well known fact for similar problems that 
increasing the degree of polynomials involved in the approximation (p-extension) is more efficient 
than refining the mesh (h-extension), see fH [3J/. 


4.1 Case A 

In Figure[l]we present the first two eigenfunctions for 7 = 1,10, 50,100,120 on [-R/2, Rj 2], We can 
observe the localization of the eigenfunction and the exponential decay far away from the wells of 
the potential. Looking at the symmetry of the eigenfunction, we see that the computation becomes 
problematic for 7 > 100 (the eigenfunctions for 7 = 120 are neither symmetric nor antisymmetric). 

I 11 Figure [2j we plot the convergence of the positive eigenvalues: 7 1 —> Aj( 7 ) and log 7 1 —> 
log 10 Aj( 7 ). We clearly observe the tunnelling effect, with the first eigenvalue Ai(y) converging 
exponentially fast to 0 while the higher eigenvalues remain of order 1 . 
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7 = 1 10 50 100 120 





Figure 2 : Convergence of the lowest eigenvalues with respect to 7 . 



Figure 3: The diffusion and drift coefficients as functions of 7 . 
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In Figure [3] we present the behavior of the diffusion and drift coefficients as functions of 7 . 
The first two columns give the result of our algorithm Step l.-Step 5. In this last column, we 
compare the numerical drift coefficient with the value /C* given by formula (4.3). More precisely, we 
display the relative error as a function of 7 . It validates the accuracy of the algorithm. 

As expected from the results on the exponential decay of the second eigenvalue, the existence of a 
spectral gap and the formulas for T> and /C, the diffusion and the drift coefficients grow exponentially 
fast as 7 —>• 00 . 









(b) 7 = 10 




40 50 

(c) 7 = 50 




Figure 4: Convergence as a function of N : D(N ), log 10 , k(N), log 10 vs. N. 


Figure [4] presents the convergence of the algorithm with respect to the number of eigenmodes 
for several values of 7 : 7 = 1 ,10, 50. Using relations (4.1), we represent 


N !->• D(N), N i->- log 10 - ——-—N i->- k(N), N i->- log 10 ^ ^ 


V * 


/c* 
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where we use the shorthand notation D(N) = for fixed R and fi (respectively k(N) = k r, ^ ,n ) 

and "D* is a reference value for the diffusion coefficient obtained for N = 50 and /C* is as above 
a numerical approximation of (4.3) with a composite rectangular rule. We observe that only 10 
eigenmodes are enough to calculate accurately the diffusion coefficient and 15 eigenmodes for the 
drift coefficient. We can conclude that our numerical method leads to an efficient and accurate 
calculation of the drift and diffusion coefficients. 


4.2 Case B — nonsmooth potential 

For the Case B, we recall that 


W{v) = ^v A -\ |u | 3 with 7 >0. 

We show similar computations as presented in the previous subsection. Figure [5] shows that the 
computations are no longer accurate beyond 7 > 7: the numerical eigenfunctions have lost their 
symmetry. For the nonsmooth potential, we can not take a larger value of 7 . 

7 = 1 7 = 5 7 = 7 





Figure 5: First two eigenfunctions on [-R/2, R/2] for 7 = 1,5,7. 




Figure 6 : Convergence of the eigenvalues. 

I 11 Figure [ 6 ] we present the convergence of the positive eigenvalues. As for the previous potential, 
we observe the exponential decay of Ai(y) to 0 as 7 increases. For large values of 7 we note a change 
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7 logi 0 X>(7) 


7 logm/C ( 7 ) 


|x:.(7)-x:(7)I 

M 7 ) 



Figure 7: Behavior of the diffusion and drift coefficients as a function of 7. 
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(a) 7 = 1 






Figure 8: Convergence as a function of N : D(N ), log 10 ^ , «(iV), log 10 vs. At. 
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of slope in the log-log graph and the loss of monotonicity for the eigenvalues A 3 , A 4 , A 5 . A further 
investigation would be necessary to decide whether this behavior is a numerical artifact or whether 
it is the actual behavior of the eigenvalues. 

I 11 Figure [7j we display the behavior of the diffusion and drift coefficients as functions of 7 
showing the exponential growth as 7 —> 00 . The comparison with (4.3) justifies the quality of the 


approximation. In Figure [ 8 ] we present the convergence of the drift and diffusion coefficients with 
respect to the number of modes for 7 = 1 and 7 = 5: as in Case A, the coefficients are well captured 
with a few eigenmodes. 


4.3 Case C — smooth potential with a linear drift 

Now we present our computations for the nonsymmetric potential 

W(v) = — l v 2 — Sv with 7 , 5 > 0 . 

The compatibility condition, necessary so that we can apply Fredholm’s alternative, does not hold: 
here 

V := f vM(v) du 7 ^ 0. 

Jr 

However, we can adapt the asymptotics in order to handle this situation where the flux of the 


equilibrium state does not vanish. Starting from (2.6), we set 


f £ (t,r,v) = f e ( t,r + ^j=,v ) , 


with f £ solution of (2.6). We check that 


d t f + ((u - V)V r / £ + div„((V r tf ★jS*)/ 6 )) = -md, 

with p £ = f f £ dv. We perform the Hilbert expansion on f £ as in Section 2 J we still have Q(/ (0) ) = 0 
and thus = p(t,r)M(v ) at leading order, while is defined by inverting 

W (1) ) = (v ~ V)V r /<°> - div„((V r tf *p)7<°>) 

= (v - V)MV r p - M ■ (V r H*p)p. 

u 


Therefore, Equations (2.8) become 


Q(xi) = (v - V)iM(v ) , 

1 dW 

Q(A) = #S W« (»)■ 
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The compatibility condition is now satisfied and we get / (1) = yy V r p + , 0- (V r f/*p)p. Considering 
the e° terms in the expansion, we finally arrive at 


d t p — V r • (T>V r p + lC(V r U -k p)p) = 0 , 


where the diffusion and drift matrices are given by the following analog of ( 2 . 11 ): 


V = 


(v — V) < 8 ) X du and /C = 


(v — V) <g) ip ch? . 


The analysis of this case can be performed as in Appendix [Aj see |23j for a similar problem. In 


5=1 5 = 5 5= 10 



(a) 7 = 1 


5=1 


5 = 5 


5 = 10 




(b) 7 = 10 



Figure 9: First two eigenfunctions on [-R/2, R/2], 


Figure [9] we present the first two eigenfunctions for 7 = 1,10 and 5 = 1, 5,10. We observe that the 
eigenfunctions are no longer symmetric. In Figure [To] we plot the first positive eigenvalue Ai(y, 5) as 
a function of 5 for several values of 7: 7 = 1 , 5, 10 , 25. As in the previous cases, Figure [XT] gives the 
approximation of the diffusion and drift coefficients using our algorithm and the comparison with 
(4.3). We observe, even for this nonsymmetric case, the good approximation of the drift coefficient. 
Figure [TX] illustrates the convergence of the algorithm for 7 = 1 and <5=1,5, 10 . As expected, we 
need more eigenmodes than for the smooth symmetric potentials in order to obtain an accurate 
approximation of the drift and diffusion coefficients. 


22 






















14 


12 

10 



0 2 4 6 8 10 


Figure 10: Ai( 7 , 5) vs. 5 for 7 = 1,5,10, 25. 





Figure 11: Diffusion and drift coefficients as functions of S for the nonsymmetric potential (??), 
with 7 = 1,5, 10, 25 (with legend as in Figure [To]) . 


5 Conclusions 

A new method for calculating the drift and diffusion coefficients in the diffusion approximation 
for a swarming model was presented in this paper. Our method is based on the calculation of the 
eigenvalues and eigenfunctions of an appropriate Schrodinger operator. This operator was obtained 
after a unitary transformation of the Markov generator that appears in the Poisson equation aris¬ 
ing in the definition of the coefficients of the limiting problem. The eigenvalue problem for this 
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(a) 7 = 1 , S = 1 



(b) 7 = 1, S = 5 



(c) 7 = 1, S = 10 


Figure 12: 
l0gl 


3lO 




Convergence as a function of N for the nonsymmetric potential (??) 
. k(A'), log 10 VS. N. 


D(N), 


Schrodinger operator was solved using a high order Finite Elements approximation. Our numerical 
method was tested to a few simple potentials and the effects of a tilt and of lack of smoothness 
of the potential on the drift and diffusion coefficients were investigated. We also investigated the 
difficulties related to the tunnelling effect that appears in the ’’semiclassical” limit 7 > 1 . 

We believe that the numerical method developed in this paper can be applied to the calculation 
of effective coefficients in a wide variety of diffusion approximations, coarse-grained and mean held 
models that appear in kinetic theory or in homogenization theory. The crucial observation is that 
in many different settings effective coefficients are given in terms of the solution to an appropriate 
linear Poisson equation pQ EH 132103102]. Thus, we believe that the spectral approach advocated 
in this paper can be of more general interest. Moreover, this approach also opens to relevant 
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perspectives: 

1. It would be interesting to develop a detailed analysis of rates of convergence, and a careful 
study of the computational cost, depending on the numerical parameters of the approximation 
(R, /i, etc). Depending on the underlying operators, the method might also benefit from the 
use of appropriate preconditioners. 

2. It would be interesting to compare the performance of the proposed numerical method with 
alternative techniques such as Monte Carlo methods 05J expansions in orthogonal polynomi¬ 
als [46] etc. 

3. The method can be extended to more general types of Poisson equations (and the calculation 
of the corresponding effective coefficients) including hypoelliptic operators of Schrodinger type 
that appear in, e.g. [26) 27], or auxilliary equations, possibly depending on the space variable, 
that appear in homogenization theory mm- 

A Analysis of the diffusion asymptotics 

A.l Set up of the problem 

In this Section we provide a few hints about the analysis of the asymptotic regime e —> 0 of the 
problem 

d t f + ^{v- V r f e - V r d> £ • V„n = -Qf , 

\j£ £ 

where on the one hand 

$ £ (t,r) = U *p £ (t,r), p £ {t,r ) = J f(t,r,v)dv, 

and on the other hand 

Qf = V v ■ (V v Wf + V J). 

For the sake of simplicity, here and below, the scaling parameter 9 is set to 1. The system is 
completed with the initial data 

f £ \ = f £ 

J 11=0 ^ * 

We set up the technical assumptions on the potentials W and U as follows. As we mentioned 
earlier, a typical example motivated from the modelling of swarming leads to the standard quartic 
potential with V v W(v) = u(^|r?| 2 — 1). Hence, in what follows, we assume that W is smooth, 
radially symmetric and coercive: 

v 1 —v W{y) depends only on |u|, is of class C 2 , and there exist C, R, a > 0 such that 
W(v) > C( 1 + |u| 2+q ) for any |u| > R, 

\V v W(v)\ < C(1 + W( v)) for any v G R d . 
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Consequently, W is bounded from below; possibly at the cost of adding to W a constant, we assume 
that 

W{v) > 0, for any v G R d . 

For U we suppose 

V r U G L°°(R d ). 

ft is likely that this assumption is far from optimal, but the technicalities to relax it are beyond 
the scope of this paper. It is convenient to set 

M(v) = Z- X e- W( ~ v) G L x (R d ), Z = j e~ w{v) dv. 

We further assume that 

W[v)M[y) G L l (R d ). 

Note that, due to the fact that W is radially symmetric, 


vM(v) dn = 0 


(A.l) 


holds. This is crucial to the analysis. As mentioned above, the discussion is strongly inspired by the 
study of the Vlasov-Poisson-Fokker-Planck system HHE9]. In this Appendix we present a proof of 


Theorem 2.1 Naturally, the proof of this theorem consists in two parts: the derivation of a priori 
estimates presented in Section [A. 2 and the passage to the limit, presented in Section A.3 


A.2 A priori estimates 

We start by observing that 


— j / f £ dv dr = 0 


holds, which means that the total mass is conserved. Here and below, we always assume 

SUPlI/lnitlU 1 < oo. 

£>0 

Owing to the regularity of U, we observe that 


verifies 


V r <F £ (t, r) = j V r U (r — r') f £ (t,r',v) dv^J dr' 


|V r $ £ (i,r)| < llVrE/IUI/^JUi. 


(A.2) 


26 







Next, we compute 


— II (W + ln(/ e ))/ e dr dr = —— // |' V v Wy/f^+2V v ^/J £ \ 2 dr dr j, // V r d> £ f £ -V v W dv dr. 


The last term can be rewritten as 


— / / V r $7 e ■ V v W dv dr = 


V r 4> 


/ fe 


V v WyfT £ + 2V,v / 7 ? 
Vi 


dr; dr 


which is dominated by 


" ,2 II fdvdr+—ll V,.W,/p - 2V„,//' =rt«dr 

< 7,rO/Uli. 7 // |v„w'VF + 2v„VF| 2 d 1 .dr, 


bearing in mind (A.2). 

In order to control the behavior at infinity of the particle distribution function, we also need to 
evaluate 

[[ \ r \f £ dr = —jn [[ vf 6 • T^-j- dr dr 


dt 




f £ — p £ M r 

v -^-• —r dr dr 

£ r 


+ pjM) ■ T dr dr 
V 1 - rl 


< 


< 


IVT^+V^rdrdr) 7 f ^? MLdvdr] 


1/2 




2 // ^ ^ + ffM ^ dU dr ) ' (/ ^ ^ dr 


1/2 


We note that the dissipation term recasts as 


I V„IT \fj £ + 2V„ v / 7| 2 dr dr = - 


, f £ 
V W — 
,,v M 


M dr dr dt. 


By the logarithmic Sobolev inequality, see [36, Th. 8.14], there exists A > 0 such that 


/e >,y f‘ 


p £ M \p £ M) p £ M 


f £ 


+ 1 j p £ M dr < A 


, f £ 
V \ — 
M 


M dr. 
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Additionally, we shall make use of the elementary inequality 


|\fs — 1| 2 < s ln(s) — s + 1 
which holds for any s > 0. It follows that 


p £ M 

Hence, we obtain 
d 


f £ 12 

■' - 1 


>Mdr = f | dv < ^ I \V v W^/ r F + 2V v y/T e \ 2 dv. 


— jj \r\f £ dv dr < JJ v 2 f £ dvdr + JJ p £ Mdrdr)+ — Jj \V v W\/J £ + 2V v \/J £ \ 2 dv dr. 

However, the assumption on the potential W guarantees that we can hnd 0 < S < oo such that 


v 2 M dv < H, 


v 2 f £ dr dr < H / (1 + W)f £ dr dr. 


In what follows, we will still denote by S a positive constant which only depends on the data, but 
which is uniform with respect to £ and t. Finally, we arrive at 


dt 


/ e (l+W+|r|+ln(/ e )) dr dr+ — / / \V v W^/J £ +2'V v \/J £ \ 2 dv dr < e(j.+ / / W/ £ drdr) 


We now use the classical trick 

s\ ln(s)| = sln(s) — 2s ln(s)l 0 < s < e -n — 2sln(s)l e -n <s < 1 < sln(s) + + 2fis, 

with s = f £ (t,x,v), H = |(|r| + W(y)). We deduce the following estimate 
W + \r\ .\ . 1 rt 


f £ ( 1 + ^ n + \Hn\)d V dr + - // \V v Wy/T e + 2V v y/p\ 2 dvdrds 


< 


\r\ + W(v) 


f £ (j. + W + |r| + ln(/ e )) dr dr 

+ J J jj \V v Wy/J^+ 2V v y/J^\ 2 dvdrds + 'd. JJ exp 
< JJ /init (l + W + |r| + ln(/f nit )) dr dr 

+z(t + JJ ex p ( - dr dr) + E J JJ Wfdvdrds. 

It remains to appeal to the Gronwall lemma to conclude with the following statement. 


dr dr 



Proposition A.l We assume that /f nit > 0 satisfies 

sup [[ (1 + |r| + W + | ln(/f nit D) /f nit dv dr = M 0 < 0. 

£>0 J J 

Let 0 < T < oo. TTien, f/iere exists 0 < Ct < oo, depending only on W, M 0 and T such that 

sup sup f f (l + |r| + W + | ln(/ e )|)/ e dr dr < C T , 

£>0 o <t<T J J 


sup / 

£>0 Jo 


12V,v / 7 F + 


dv dx d t < Ct- 


This estimate can be translated by means of compactness properties, as an application of the 
Dunford-Pettis theorem, see e. g. [22., Section 7.3.2 & 7.3.4], 

Corollary A.l The sequence (f £ ) e>0 is weakly com-pact in L 1 ((0,T) xl^x R d ). 

Corollary A.2 We can write f £ (t,r,v ) = p £ (t,r)M(v) + y/eg £ (t,r,v ) ; with f g £ dv = 0, where the 
sequence (p £ ) e>0 is weakly compact in Id^O, T) x R d ) and the sequence (g £ ) £>0 is weakly compact 
in L 1 ((0,T) xR d x R d ). 


Proof. The compactness of the macroscopic density p £ is a direct consequence of Corollary A.l 
Next, let A C (0, T) x R d x R d be a measurable subset. Reproducing a manipulation already 
detailed above we get 



\g £ \dvdrdt = III drdr dt 


A 


< 2 




IA 

\ 1/2 

(f £ + p £ M ) dr dr dt) 

<[ 2 I I j\f £ + p £ M) dr dr dt^ 

< -Jwi (JJJ if' + (CM) dt> dr d t) 1 


\AEtTZE 1 AvirAt C 


£ J 

\2V v VF + V v W^f £ \ 2 


dr dr dt 


1/2 


Coming back to Corollary A.l it proves the equi-integrability of the sequence g £ . 


Owing to the behavior of W for large |r|’s, we shall use the fact that the reasoning can be 
applied to f £, f>(v) with tests functions r H» if(y) verifying lim^i^oo = 0 and to g £ if( r) as well, 
with tests functions r H- ip(v) verifying lim|„|_ 5 . 00 = 0 In particular, using (|A. lh , we have 

\/W(v) - 
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Corollary A.3 The sequence defined by 


J £ (t,r) 


f £ (t,r,v)dv 

sje 


J v g £ [t , r, v) dr; 


is weakly compact in L 1 ((0,T) x M d ). 


A.3 Passage to the limit 

According to the discussion above, we can assume, possibly at the price of extracting a subsequence, 
that 

f £ f weakly in ^((0, T)xK d x R d ). 

It means that 


lim 

£— 



f £ (t, r, r, v ) dr dr df 



f(t, r, r, v ) dr; dr dt 


holds for any trial function G L°°((0,T) xW d x R d ). In fact, owing to the uniform estimate on 
the integral of Wf £ the convergence applies for fi(t,r,v) = ((t, r)0(r;), with ( G L°°((0,T) x R d ) 
and </> G C 0 (M d ) such that lim^i^oo = 0. Accordingly, we have 


P 


£ 


P 


J fdv weakly in L 1 ((0,T) x R d ). 


Integrating the equation with respect to v, we obtain 


d t p £ + V r ■ J £ = 0. 


(A.3) 


With Corollary A.2 and Corollary A.3 we can also suppose that 


g £ g weakly in L 1 ((0,T) xR d x R d ), 
J e —*■ J weakly in L 1 ((0, T) x R d ). 


Letting e go to 0 in (A.3) we are thus led to 

dtp -\~ V r • J — 0. 


Furthermore, for any trial function Q G VF 1,0O (R d ), (A.3) and Corollary A.3 imply that 


p £ {t + h, r)((r) dr — / p £ [t, r)C(r) dr = 


pt-\-h 


J £ (s,r ) • V r C(r) drds 
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can be made arbitrarily small as h goes to 0, uniformly with respect to e. Up to a suitable regu¬ 
larization argument, we deduce that p £ is compact in C°([0, T]; L 1 (M d ) — weak). As a consequence, 
extracting further subsequences if necessary, we deduce that 


V r <3> £ (t, r) = / V,.f/(r — r')p £ (t, r') dr'—G-VV^t, r) = / V r Z/(r — r')p(t,r') dr' 


holds a.e. (0, T) x M , with the uniform estimate (A.2). 


It remains to identify the limiting current J. Going back to the proof of Corollary |A.2 
justify that 


we 


lim 

£—>-0 


g £ (t, r, r)-0(t, r, v) dr dr dt = 


g(t, r, r)0(t, Ar) dr dr dt 


holds for any trial function -0 G L°°((0,T) x x M d ) as well as for 0(t, r,r) = £(t,r)0(r), with 
C G L°°((0, T) x M d ) and 0 G C'°(M d ) such that lim^i-Hx) —^== = 0. In particular, we have 


J e (t,r) = / rc£(t, r, r) dr — 1 J(t,r ) = / vg(t,r, r) dr weakly in L 1 ((0,T) x M d ). 


Then, for any G C“((0, T) x M“) and 0 G we have 


f E (t,r,v ) (f>(v)d t ((t,r) dr drdt 


fT 


/ £ (t, r, r)0(r)r ■ \ / r ((t,r) dr dr dt 


/ £ (t,r, r) V„0(r) • V r <l? £ (t, r) £(t,r) dr drdt 


——j= / £0(r) / £ (t, r, r)C(t,r) dr dr dt = / £0(r) c/ £ (t, r, r)£(t, r) dr dr dt, 


'0 


'0 


where £ is the Zr-adjoint operator to Q defined in (2.5). Letting e go to 0, we are led to 

rT 


'0 


£0(r) g(t, r, r)C(t, r ) dr dr dt = 


'0 


+ 


/(t, r, r) 0(r)r ■ V r 0(t, r) dr dr dt 


/(t, r, r) V„0(r) ■ V r <h(t, r) £(t,r) dr drdt. 


For the nonlinear term, we have combined the weak convergence of f £ in L 1 and the pointwise 
convergence of the uniformly bounded sequence V r <F £ , see |22j Lemma 7.62], By virtue of (A.2), 
we have 

/ (t, r, r) = p(t,r)M( r). 
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Therefore, g is characterized by the relation 


rT r r rT 

£<p(v) g(t, r, v)((t, r ) dv dr dt = — 


4>{v)vM{v) dr ) ■ p(t,r)V r C(t,r) dr dt 


V V (p(v) M (v) dv ) ■ V r $(f,r)p(f,r)C(f, r)drdt, (A.4) 


for any ( G C£°((0, T) x M d ), </> G C£°(R d ). Note in particular that 

[ V v (f>{v) M(v)dv= [ <f>(v)V v W{v)M(v) dv. 


Going back to (A.4), it follows that 

r T 


C</)(v)g(t,r,v)((t,r)dvdrdt < \\(f)\\L\Mdv) WpWl^^l^r^T 
x(||V r C|Uao||u>/M|| L2 + ||V r Z/|Ucc||p|| Loo(0iT \\(\\ L °o\\V v Wy/M\\ L ^ 

which makes the regularity of Qg precise. 


We recall that Q is a differential operator associated to a natural quadratic form with domain 

and we have 

-jQfg dv = f Vu ~~ ■ Mdv. 

Note that L 2 (-^-r) C L 1 , and by the Sobolev inequality [3J Corollary 2.18], we check that Q is 
coercive on the closed subset of functions with zero-mean. Therefore, the Fredholm alternative 
applies as follows. 


Lemma A.l For any h G y) such that f hdv = 0 (resp. h G L 2 (M dr>) such that f hMdv = 

0), there exists a unique f G D(Q) such that Qf = h and f f dv = 0 (resp. f G D(C) such that 
Cf = h and f fMdv — 0). 


Consequently, we can define x solution of C~x 
rewrite the limiting current as follows 


G L 2 (M dv) as in (2.12). Therefore, we can 


J £ (t,r ) 


J vg £ (t,r,v) d v J(t,r) 


j vg(t,r,v ) dv 


£x(v) g(t,r,v) dv. 
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Identifying limits with (A.4), np to a regularization argument, we arrive at 


x(v) ® vM{y) dn ) pit, r) V r C(£, r) dr dt 


x(v) <g> V v W(v)M{y) dv I V r $(t, r)p{t, r ) dr dt. 


I J J(t,r)((t,r)drdt = -j 

r T r 


It corresponds to the dual formulation of the expected relation. Indeed, we observe that, on the 
one hand 

/ WO ® vM{v) dv = Ji(v)®Qx(v)dv = Jv® X (v)dv = -V, 
and, on the other hand 

/ m 9 V„W(v)M(v) dv = / m ® <M«0 to = /»® «<«o * = -«■ 

Hence we have obtained 


J = — ( W r p + /CV. r $ p). 


It ends the proof of Theorem 2.1 


While this is not necessary for establishing the connection between the kinetic model and the 
drift-diffusion equation in the regime e —y 0, it is possible to improve the compactness of the 
macroscopic density from weak to strong. The proof relies on the combination of a renormalization 
argument and velocity averaging techniques m This is detailed in mm • We start with the 
following averaging lemma. 

Lemma A.2 UUA Prop. 4-1] Let h £ be a uniformly bounded sequence in L 2 ((0,T) x x M d ). We 
assume that 

{V~ed t + v ■ X7 r )h £ = h £ 0 + V„ ■ h\, 

with h £ 0 and h\ bounded in T 1 ((0, T) x M. d x M d ). Then, for any function if 6 we have 


sup 

£>0 

Let 0 < 5 < 1. We set 


(, h £ {t , r + h, v) — h £ (t, r, u)) if(v) dn 


Ps(z) = 


L l I^H*'0 


-> 0 . 


1 + Sz 

Therefore h £ = fds(f £ ) is bounded in L 1 flL oo ((0, T) xM d xM d ), thus in L 2 ((0, T) xM d xM d ), uniformly 


with respect to e (but the bound depends on 5). Hence Lemma A.2 applies to h £ = (3s(f £ ) for any 
0 < S < 1 fixed with 

* = punr (- + E . 
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and 


We deduce that 


hi = v r * e ftcn-^cnv^ 


V v Wy/p + 2V. uV Ap 


sup 

£>0 


(Ps{f e ){t,r + h,v) - p s (f)(t,r,v)) ip(v) dv 


-» 0 

Li PI— 


holds for any 5 > 0. This property passes to f £ owing to the equi-integrability of the sequence. 
Indeed, we split 

J {F(t,r,v) - p s (f £ )(t,r,v)) ^(v)dvdrdt 

by considering separately z = ( t,r,v ) such that \z\ > M, or \f £ {z)\ > A and the complementary 
sets. The integral of the former can be made arbitrarily small by chosing M, A large enough. The 
integral of the latter is dealt with by using the fact that /3s(s) converges uniformly on compact sets 
to s. We arrive at 


sup 

£>0 


(/ £ (f, r + h, v ) - f £ (t, r, v)) dv 


->■ 0. 

jji 1^1 — 


Note that owing to the weighted estimates on f £ we can remove the restriction of considering 
compactly supported trial functions. This is finally combined to (A.3): the a priori estimates tells 
us that d t p £ is bounded in L 1 (0, T; W _1,1 (M d )). We conclude by using a standard approximation 
argument that p £ satisfies the Weil-Kolmogorov-Frechet criterion, see [22, Theorem 7.56]. Therefore 
p £ converges to p strongly in T 1 ((0, T ) x M d ). 
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